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ABSTRACT 

The isothermal gravitational collapse and fragmentation of a region within a molec- 
ular cloud and the subsequent formation of a protostellar cluster is investigated numer- 
ically. The clump mass spectrum which forms during the fragmentation phase can be 
well approximated by a power law distribution dN/dM oc M~ L5 . In contrast, the mass 
spectrum of protostellar cores that form in the centers of Jeans-unstable clumps and 
evolve through accretion and iV-body interactions is described by a log-normal distri- 
bution with a width that is in excellent agreement with observations of multiple stellar 
systems. 

Subject headings: ISM:kinematics and dynamics — ISM:structure — methods:numerical 
— stars:formation 
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1. Introduction 

Understanding the processes leading to the 
formation of stars is one of the fundamental chal- 
lenges in astronomy and astrophysics. With the 
advent of new observational techniques and in- 
struments, especially in the IR and radio wave- 
bands, a vast amount of astronomical data about 
star forming regions has been accumulated. How- 
ever, on the theoretical side not much progress 
has yet been made. Analytical models of the 
star formation process are restricted to describ- 
ing the collapse of isolated, idealized objects (for 
an overview see Whitworth & Summers 1985). 
Much the same applies to numerical studies (e.g. 
Bonnell & Bastien 1993, Boss 1997, Burkert et 
al. 1996, 1997, Nakajima & Hanawa 1996). Star 
formation is a complex self-gravitating, magneto- 
hydrodynamical problem, which includes the ef- 
fects of heating and cooling, and feedback pro- 
cesses from newly formed stars. Furthermore, it 
is influenced by the galactic environment. Taking 
into account all these processes with high spatial 
resolution exceeds by far present computational 
capabilities. 

Previous numerical simulations of the collapse 
and fragmentation of molecular cloud regions 
have shown that a large number of condensed ob- 
jects can indeed form on a dynamical timescale as 
a result of gravitational fragmentation (e.g. Lar- 
son 1978, Monaghan k, Lattanzio 1991, Keto 
et al. 1991). In these studies, the clouds were 
treated as isolated gaseous spheres which col- 
lapsed completely onto themselves. Instead, we 
study a small region embedded in a large, stable 
molecular cloud complex where only the over- 
dense regions are able to contract due to self 
gravity. We assume the molecular cloud is sup- 
ported on large scales by turbulence and/or other 
processes. Previous numerical models were also 
strongly constrained by numerical resolution. Lar- 
son (1978), for example, used just 150 particles in 
an SPH-like simulation. Whitworth et al. (1995) 
and Turner et al. (1995) were the first who ad- 



dressed star formation on larger scales in detail 
using high-resolution numerical models. How- 
ever, they studied a different problem: fragmen- 
tation and star formation in the shocked inter- 
face of colliding molecular clumps. While clump- 
clump interactions are expected to be abundant 
in molecular clouds, the rapid formation of a 
whole star cluster requires gravitational collapse 
on a larger scale which contains many clumps 
and dense filaments. 

In this letter, we extend previous studies of 
the collapse of isolated objects to the regime of 
the isothermal collapse and fragmentation of a 
gravitationally unstable region embedded in the 
interior of a molecular cloud. We present a high- 
resolution numerical model of the dynamical evo- 
lution and follow the fragmentation into dense 
protostellar cores. The temperature and the den- 
sity is chosen such that the region is highly grav- 
itationally unstable and forms a hierarchically- 
structured protostellar cluster. The results of 
this study, i.e. the properties of the dense clumps 
and of the newly formed protostellar cores are 
compared with observations. 

2. Numerical Technique 

To follow the time evolution of the system, 
we use smoothed particle hydrodynamics (SPH: 
for a review see Monaghan 1992) which is in- 
trinsically Lagrangian and can resolve very high 
density contrasts. The code is based on a version 
originally developed by Benz (1990). We adopt a 
standard description of artificial viscosity (Mon- 
aghan & Gingold 1983) with the parameters a v = 
1 and (3 V = 2. The smoothing lengths are variable 
in space and time such that the number of neigh- 
bors for each particle remains at approximately 
fifty. The system is integrated in time using a sec- 
ond order Runge-Kutta-Fehlberg scheme, allow- 
ing individual timesteps for each particle (Bate 
et al. 1995). Once a highly-condensed object has 
formed in the center of a collapsing cloud frag- 
ment and has passed beyond a certain density, 
we substitute it by a 'sink' particle which then 
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continues to accrete material from its infalling 
gaseous envelope (Bate et al. 1995). By doing so 
we prevent the code time stepping from becom- 
ing prohibitively small. This procedure implies 
that we cannot describe the evolution of gas in- 
side such a sink particle. However, at some stage 
of the gravitational collapse the SPH resolution 
limit (Bate & Burkert 1997) would be reached 
in the fragment anyway. For a detailed descrip- 
tion of the physical processes inside a protostellar 
core, i.e. its further collapse and fragmentation, 
a new simulation just concentrating on this sin- 
gle object with the appropriate initial conditions 
taken from the larger scale simulation would be 
necessary (Burkert at al. 1998). 

To achieve high computational speed, we have 
combined SPH with the special purpose hard- 
ware device GRAPE (Sugimoto et al. 1990, Ebi- 
suzaki et al. 1993), following the implementa- 
tion described by Umemura et al. (1993) and in 
greater detail by Steinmetz (1996). Since we wish 
to describe a region in the interior of a globally- 
stable molecular cloud we have to prevent global 
collapse. Therefore, we use periodic boundaries, 
applying the Ewald (1921) method in an PM-like 
scheme (Klessen 1997). In this letter we present 
a simulation with 500 000 SPH particles. 

3. Initial Conditions for Cloud Fragmen- 
tation 

The structure of molecular clouds is very com- 
plex, consisting of a hierarchy of clumps and 
filaments on all scales (for a review see Blitz 

1993) . Many attempts have been made to iden- 
tify the clump structure and derive its prop- 
erties (Stutzki &: Glisten 1990, Williams et al. 

1994) . We choose as starting conditions den- 
sity fields with Gaussian random fluctuations 
that follow a power spectrum P{k) oc 1/k 2 , i.e. 
large scale fluctuations have on average large am- 
plitudes whereas the amplitudes of short wave- 
length modes decay quadratically. The fields are 
generated by applying the ZePdovich (1970) ap- 
proximation to an originally homogeneous gas 



distribution: we compute a hypothetical field of 
density fluctuations in Fourier space and solve 
Poisson's equation to obtain the corresponding 
self-consistent velocity field. These velocities are 
then used to advance the particles in one big step 
5t. 

4. A Case Study 

As study, we present the time evolution 

of a region in the interior of a molecular cloud 
containing a total mass of 222 Jeans masses, de- 
termined from the temperature and mean density 
of the gas. Figure |l| depicts snapshots of the sys- 
tem initially, and when 10, 30 and 60 per cent 
of the gas mass has been accreted onto the pro- 
tostellar cores. Note that this cube has to be 
seen periodically replicated in all directions. Ini- 
tially, pressure smears out small-scale features, 
whereas large-scale fluctuations start to collapse 
on themselves, and into filaments and knots. Af- 
ter t ~ 0.9Q, the first highly-condensed cores 
form in the centers of the most massive and dens- 
est Jeans-unstable gas clumps and are replaced 
by sink particles. Soon, clumps of lower initial 
mass and density follow, altogether creating a 
hierarchically-structured cluster of accreting pro- 
tostellar cores. 

4.1. Scaling Properties 

The gas is isothermal. Hence, the calculations 
are scale free, depending only on one parameter: 
the dimensionless temperature T = Ei nt /\E pot \, 
which is defined as the ratio between the inter- 
nal0 and gravitational energy of the gas. The 
model can thus be applied to star-forming regions 
with different physical properties. In the case of a 

lr Time is measured in dimension- less units with t = being 
denned by the start of the SPH simulation. For adequate 
timing, the Zel'dovich shift interval St — 1.5 has to be 
added. As reference, the free-fall time of the isolated cube 
would be Tff = 1.4. 

2 In the case of isotropic turbulence, the non-thermal (tur- 
bulent) contributions can also be accounted for in this 
expression, E int = E^ erm + E tU rb- 
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Fig. 1. — Time evolution and fragmentation of a region of 222 Jeans masses in the interior of a molecular 
cloud with initial Gaussian density fluctuations with power law P(k) oc l/k 2 . Collapse sets in and soon 
forms a cluster of highly-condensed cores, which continue to accrete from the surrounding gas reservoir. 
At t = 1.6 about 10% of all the gas mass is converted into " protostellar" cores (denoted by black dots). 
At t = 2.0 and t = 2.8 these values are 30% and 60%, respectively. The initial number of particles used 
for the SPH simulation is 500000. For legibility only every tenth particle is plotted. 



dark cloud with mean density n(H2) ~ 100 cm" 3 
and a temperature T ~ 10 K like Taurus- Auriga, 
the computation corresponds to a cube of length 
10 pc and a total mass of 6 300M Q (with the 
Jeans mass being Mj = 28 M ). The dimen- 
sionless time unit corresponds to 2.2 x 10 6 yrs. 
For a high-density star-forming region like Orion 
with n(H2) ~ 10 5 cm" 3 these values scale to 
0.32 pc and 200 M Q , respectively, again assum- 
10 K. The Jeans mass is 0.9 Mq and the 
time scale is 6.9 x 10 4 yrs. 



4.2. The Importance of Dynamical Inter- 
action and Competitive Accretion 

The location and time at which protostel- 
lar cores form, is determined by the dynami- 
cal evolution of their parental gas clouds. Be- 
sides collapsing individually, clumps stream to- 
wards a common center of attraction where they 
merge with each other or undergo further frag- 
mentation. The formation of dense cores in the 
centers of clumps depends strongly on the rela- 
tion between the timescales for individual col- 
lapse, merging and sub-fragmentation. Individ- 
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Fig. 2. — Mass growth history for 14 selected 
protostellar cores selected from a total of 55. The 
numbers reflect the order of their formation. 

ual clumps may become Jeans unstable and start 
to collapse to form a condensed core in their cen- 
ters. While collapsing, these clumps also interact 
with each other. When clumps merge, the larger 
new clump continues to collapse, but contains 
now a multiple system of cores in its center. Now 
sharing a common environment, these cores com- 
pete for the limited reservoir of gas in their sur- 
rounding (see e.g. Price & Podsiadlowski 1995, 
Bonnell et al. 1997). Furthermore, the "proto- 
stellar" cores interact gravitationally with each 
other. As in dense stellar clusters, close encoun- 
ters lead to the formation of unstable triple or 
higher order systems and alter the orbital pa- 
rameters of the cluster members. As a result, a 
considerable fraction of "protostellar" cores get 
expelled from their parental clump. Suddenly 
bereft of the massive gas inflow from their col- 
lapsing surrounding, they effectively stop accret- 
ing and their final mass is determined. Ejected 
objects can travel quite far and resemble the 
weak line T Tauri stars found via X-ray obser- 
vation in the vicinities of star-forming molecular 
clouds (e.g. Neuhauser et al. 1995, Wichmann et 
al. 1997). 

In Fig. [2], we plot the accretion history for 



14 representative protostellar cores in the sim- 
ulation. The objects are numbered according to 
their time of formation. The figure illustrates 
the following trends: (a) The cores which form 
first tend to have the largest final masses. They 
emerge from the initial clumps with the high- 
est densities. The collapse of these large ini- 
tial clumps contributes a considerable fraction 
of their total mass, (b) Matter that contracts 
into dense cores at later times (say t > 2) has 
already undergone considerable dynamical evo- 
lution. Small initial clumps stream towards each 
other along filaments. At the intersections of fil- 
aments they merge and may undergo rapid col- 
lapse, when enough mass is accumulated, (c) 
Once dense cores have formed, they evolve due to 
accretion, competing for gas from the surround- 
ing reservoir and interacting dynamically, as de- 
scribed above. For example, at t pa 1.8, core #19 
is expelled from a dense clump at the intersec- 
tion of two massive filaments by a triple interac- 
tion with cores #1 and #17. It stops accreting. 
However, it still is bound to the gas knot which 
grows in mass due to continuous infall. It falls 
back onto the clump of gas and resumes accreting 
at t pa 2.0. Cores #9 and #41 are also expelled 
from their parental clumps but, unlike core #19, 
their accretion is terminated completely. These 
dynamical interactions between cores are an im- 
portant agent in shaping the mass distribution. 

4.3. Mass Spectrum — Implications for 
the IMF 

Figures ||a - d describe the mass distribution 
of identified gas clumps (thin lines) and of pro- 
tostellar cores (thick lines). To identify indi- 
vidual clumps we have developed an algorithm 
similar to the method described by Williams et 
al. (1994), but based on the framework of SPH. 
As a reference, we also plot the observed canoni- 
cal form for the clump mass spectrum, dN/dM oc 
AT -1 - 5 (Blitz 1993), which has a slope of -0.5 
when plotting N versus M. Note that our ini- 
tial condition does not exhibit a clear power law 
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Fig. 3. — a) - d) Mass distribution of gas clumps (thin lines) and of protostellar cores (thick lines) at 
times t = 0.0, 0.7, 1.3 and 2.0 when 0%, 10%, 30% and 60% of the total gas mass is condensed in cores, 
respectively. The vertical lines indicate the resolution limit of the simulation with 500 000 particles, and 
the dashed lines illustrate the observed clump mass spectrum with dN/dM oc M" L5 (Blitz 1993). e) 
Comparison of the final core mass spectrum (thick line) with different observationally based models for 
the IMF. The thick dashed line denotes the log-normal form for the IMF, uncorrected for binary stars 
as proposed by Kroupa et al. (1990). In order for the peaks of both distributions to overlap, a core star 
formation efficiency r\ has to be assumed. The agreement in width is remarkable. The multiple power-law 
IMF, corrected for binary stars (Kroupa et al. 1993) is shown by the thin solid line. As reference, the 
thin dashed line denotes the Salpeter (1955) IMF. Both are scaled to fit at the high-mass end of the 
spectrum. All masses are scaled to the overall Jeans mass in the system. 



clump spectrum, but instead consists preferen- 
tially of small scale fluctuations. However, these 
are quickly damped by pressure forces and dur- 
ing the subsequent non-linear gravitational col- 
lapse a power-law mass spectrum is formed with 
a slope that is similar to the observed clump mass 
spectrum (Fig. ||b,c). In all panels, the vertical 
lines indicate the SPH resolution limit for 500 000 
particles (Bate & Burkert 1997). 

A common feature in all our simulations is 
the broad mass spectrum of "protostellar" cores 
which peaks slightly above the overall Jeans mass 
of the system. This is somewhat surprising, 
since the initial fluctuations span a vast range 
of masses and peak densities, and the evolu- 
tion of each core is heavily influenced by com- 
plex merging and collapse processes. However, 
in a statistical sense, the system retains "knowl- 



edge" of its (initial) average properties. The 
present simulations cannot resolve subfragmen- 
tation in condensed cores. However, detailed 
simulations show that perturbed cores tend to 
break up into multiple systems (e.g. Burkert et 
al. 1996, 1997). Here, we can only determine 
the mass function of multiple systems without 
breaking them down into a mass function of sin- 
gle stars. Our simulations predict an initial mass 
function with a log-normal functional form, if 
the mass of a multiple system that forms within 
a condensed core is roughly proportional to the 
core mass, M* = 77 x M core . Figure ||e compares 
the results of our calculations with the observed 
IMF for multiple systems (Kroupa et al. 1990). 
The maximum of the observed IMF is located 
at 0.23 M . The calculated core mass distribu- 



tion peaks at ~ 2Mj. In Section 4.1, we found 
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that when we scale the simulation to the con- 
ditions in low-density (e.g. Taurus) and high- 
density star- forming regions (like Orion), we ob- 
tain Jeans masses of 28 M Q and 0.9 M , respec- 
tively. Thus, to reproduce the observed peak 
in the IMF, we take r] ~ 0.005 for Taurus and 
rj ~ 0.125 for Orion, i.e. the implied star forma- 
tion efficiency for a low-density Taurus-like re- 
gion is much lower than for cloud regions of high 
density like Orion. With these star formation 
efficiencies, the agreement between the observed 
IMF for multiple sytems (thick dashed line; val- 
ues from Kroupa et al. 1990) is excellent. Note 
that although the peak depends on our choice of 
77, the agreement in the width of the distribution 
does not depend on this scaling. For compari- 
son, the IMF corrected for binary stars (Kroupa 
et a. 1993) is indicated as thin solid line, together 
with the mass function from Salpeter (1955) as 
thin dashed line. 

5. Summary and Discussion 

Since collapse and fragmentation in molecular 
clouds is an extremely complex and dynamical 
process, many authors have sought to understand 
the stellar initial mass function as resulting from 
a sequence of statistical events which may natu- 
rally lead to a log-normal IMF (see e.g. Zinnecker 
1984, Adams k Fatuzzo 1996; also Price & Pod- 
siadlowski 1995, Murray & Lin 1996, Elemegreen 
1997). However, using numerical simulations, it 
is possible to identify underlying processes which 
may contribute to the form of the stellar ini- 
tial mass function. In the calculations presented 
here, we find several trends. The "protostel- 
lar" cores that form first are generally formed 
in the clumps with the highest initial density, 
and tend to have the highest final masses. Cores 
that form later, form from gas that was initially 
in low-density clumps or distributed gas which 
converged to form a higher-density clump be- 
fore quickly collapsing. Overlaid on these gen- 
eral trends, dynamical interactions between in- 
dividual cores can act to terminate accretion on 



to a core by ejecting it from a clump, thus set- 
ting its final mass. The excellent agreement be- 
tween the numerically-calculated mass function 
and the observed IMF for multiple stellar sys- 
tems (Kroupa et al. 1990) strongly suggests that 
these gravitational fragmentation and accretion 
processes dominate the origin of stellar masses. 
In a subsequent paper, the results from calcula- 
tions spanning a larger range of the parameter 
space relevant for molecular clouds shall be dis- 
cussed in detail. 
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